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with longitudinal expansion 



CN \ Yoshitaka Hatta a and Akihiro Nishiyama 6 

a Faculty of Pure and Applied Sciences, University ofTsukuba, Tsukuba, Ibaraki 305-8571, Japan 



> 



o 



X 



b 



Maskawa Institute for Science and Culture, Kyoto Sangyo University, Kyoto 603-8555, Japan 



-i— > 
o 

O 

^ Abstract 

Motivated by the problem of thermalization in heavy-ion collisions, we present numerical simu- 
lations of the nonequilibrium evolution of the O (N) model in 1+ 2 dimensions with longitudinal 
expansion and in the presence of a background field. We work in the NLO approximation of the 
1 /N expansion and consider both the strong and weak coupling cases. Special emphasis is put 
on the difference between the 2PI formalism and the classical statistical approach. In the 2PI 
case at strong coupling, we find some evidence of the Bose-Einstein (exponential) distribution 



in the expanding system. 



1. Introduction 

Understanding how and when thermalization is achieved in the fireball created in heavy- 
ion collision experiments is one of the most pressing, but difficult questions in physics of the 
Quark Gluon Plasma (QGP). The system starts out with a very peculiar initial condition de- 
scribed by the Color Glass Condensate (CGC) framework [1], and eventually ends up with the 
hydrodynamic regime characterized by the thermodynamic properties of QCD. The transient 
regime in between is a 'missing link' about which our knowledge is severely limited due to the 
highly nonlinear dynamics of QCD which is furthermore complicated by the one-dimensional 
expansion of the system. 

Analytical insights are available for certain aspects of the evolution, and various scenarios 
for thermalization have been put forward. [See [2] and references therein.] But if one wants 
to go beyond parametric estimates and get the whole evolution under quantitative control, nu- 
merical simulations based on the first principles seem to be necessary. So far, simulations 
based on the classical statistical approximation have proven to be very practical for this purpose 
[SAUL 0,0]. These simulations nicely demonstrate, among others, the exponential growth of 
soft gluons triggered by the decay of classical fields. However, they also indicate that in the 
weak coupling picture the unstable regime develops in an unrealistically (from the viewpoint 
of heavy-ion phenomenology) large span of time. Besides, the classical statistical approxima- 
tion is a low-momentum effective theory, so 'thermalization' achieved in this framework is the 
classical Rayleigh-Jeans type, and not the quantum Bose-Einstein distribution or the Boltz- 
mann (exponential) distribution which is often regarded as a signal of the QGP in heavy-ion 
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collisions. It is then difficult to tell whether one can form the truly equilibrated QGP on a rea- 
sonable time scale, or just some far-from-equilibrium gluonic matter which mimics the QGP 
by acquiring a constant energy-pressure relation ('equation of state') due to the phenomenon 
of prethermalization H. 

In principle, the latter difficulty can be solved by the two-particle irreducible (2PI) formalism 
which has a variety of applications in nonequilibrium physics and has been recently discussed 



in the context of heavy-ion collisions in the CGC picture [9] (see, also, llOfl). In the 2PI formal- 
ism, the Bose-Einstein distribution is guaranteed as it consistently includes the relevant terms 
responsible for the quantum equilibration. However, unlike the classical statistical approach, 
usually the collision terms in the 2PI approach are evaluated to fixed order in the coupling 
constant. While such an approximation may be justified for certain purposes, in the presence 
of instability at weak coupling it will break down in the soft sector because higher order loop 
diagrams, proportional to large powers of the particle density, are no longer suppressed at the 
height of instability. 

Given that the classical statistical approach and the 2PI formalism have their own advan- 
tages (and disadvantages), one may hope to find a unified approach which captures both these 
advantages in a complementary manner. Unfortunately such a formalism is not available in 
QCD at the moment. There is, however, one (and only one) known example where this can be 
done — It is the scalar 0(N) theory at large N [1 1]. To next-to-leading order (NLO) in the 1/N 
expansion, one can resum the 2PI loop diagrams (including the terms omitted in the classical 
statistical approach) to all orders in the coupling constant. It is thus a very attractive toy model 
from which one can learn important lessons for QCD. Numerical simulations of this model have 
been previously performed in the flat coordinates [|l2[ Hill . In this paper we report the result of 
numerical simulations of this model in (l+2)-dimensions with one-dimensional expansion (one 
time dimension, one transverse dimension, and one longitudinal (expanding) dimension) and in 
the presence of a background field. Our main focuses are the difference between the 2PI and the 
classical statistical approaches, and the impact of longitudinal expansion on the equilibration 
process. For previous 2PI simulations in (spherically) expanding systems, see Ill4ll . 

2. The model 

Our model field theory is the 0(iV)-symmetric \cp A theory in (l+2)-dimensions = 

(t,x T ,z) 

S = J dtdx T dz QfljA^O - 7^ 2 0a0a ~ ^ -{4 >a$ "a) J , (1) 

where a = 1, 2, • • • , N. A is the dimensionless coupling constant and the dimensionful pa- 
rameter cr sets the energy scale of the problem. We assume that the system is undergoing 
one-dimensional expansion in the ^-direction. It is then convenient to switch to the W-rf 
coordinates defined by 

T = Vt 2 - z 2 , r/ = tanh" 1 |. (2) 
r is the so-called proper time, and rj is the rapidity. In these coordinates, the action reads 

S = J TdTdx T dr) Q(<9 T a ) 2 - 7^(<V a ) 2 - ^(d T (f> a ) 2 - ^m 2 ^ a (f) a - ^(<f) a (f> a ) 2 j , (3) 
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and the classical equation of motion takes the form 



d 2 T + \d T -^dl-dl + m' + ^ Mb 



0. 



(4) 



We shall be interested in the evolution of the system characterized by a nonzero field expectation 
value at the initial time r = r 



(Uro)) 



>al 



6N 



A ' 



(5) 



which we take to be boost invariant (independent of rj) and, for simplicity, homogeneous (inde- 
pendent of x T ). Such a background field is to mimic the initial condition of heavy-ion collisions 
in the color glass condensate picture. In that case, weak (QCD) coupling is assumed so that the 
classical gluon fields are strong ~ 0(1/ y/a^) 3> 1. Here we relax this assumption and treat 
A as a free parameter. 

At t > t the classical field decays, and fuels energy for quantum fluctuations to grow. The 
latter is characterized by the statistical function 



F ab (x,x') = l -({4> a {x),Ux')}) - (<f> a (x))(<f> b W)) , 



and the spectral function 



p ab (x,x') = i{[<f) a (x),(f) b (x')]) . 



(6) 



(7) 



We choose the initial condition such that (<j) a (T)) = for a ^ 1 at all times, and introduce 
the 0(N — 1) -symmetric notation Fu = F\\, F 22 = F 33 = • • • = Fnn = F±. Since the 
classical field is independent of xt and 77, we work in the Fourier space and denote the conjugate 
momenta to x T and rj by p T and p v , respectively. 

The evolution of the system is described by the following set of coupled equations ('Kadanoff- 
Baym equation') for (0), F and p : 



# + \d T + m 2 + ^ ^> 2 (r) + 3F n (r, r) + g F bb (r, r) 



4>{t) 



-f 

J TO 



where 



F(r,r) = J 



dp T dpr, , v 
(2tt) 2 ^ T > PT > P ^> 



(8) 



(9) 



and (4>i (t)} is simply denoted as <f)(r). The equations for the statistical and the spectral func- 
tions are 

^ + ^ + ^+PtW + M 2 & (0) F bc (r,r',p) 

- [ T r"dr"^ ab (r,T",p)F bc (r",r',p)+ f r"dr"YF ab {ry ,p)p bc (r" ,r' ,p) , (10) 
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T T 



,2 

Pbc(r, t', p) = - I r"dr"E^(r ) t", p)p bc (r", r', p) 



d 2 T + ^d T + ^+p 2 T )5 ab + M 2 ab ( ( l ) ) 



(11) 



where the effective mass term is 

M a 2 b (0) = m 2 5 ab + ^ (F cc (r, r) + 2 ) 5 afe + ^ (F ab (r, r) + 2 5 al 5 w ) . (12) 

The self energies Ef fe and Tj P ab are evaluated in the 2PI formalism to next-to-leading order 
(NLO) approximation in the 1/N expansion. Their expressions are rather lengthy and are not 
reproduced here. We refer the reader to the original paper fUlll and, importantly for the present 



purpose, its extension [|15ll to systems with a certain class of curved coordinates including the 
t-T) coordinates. 

3. Numerical simulations: Setup 

Thanks to the manifestly causal structure, the equations (l8l-(fTTI) can be straightforwardly 
solved on a lattice. We consider a lattice of size 64 x 64 in the momentum space with the 
periodic boundary condition. The momentum and time are discretized as 

— = 7Ta = "HT> ^ = 7^ = T^T' -32 <n T ,n v < 32, (13) 

a b4a T a 64 64a,, 16 



T = na t , t = I0a t , (14) 

where we have chosen the lattice spacings a T = 2/a, a v = 0.5, and a t = OA/a. Due to the 
symmetry of the problem p <-> — p, the actual simulations may be performed on a 32 x 32 lattice 
by restricting to the region p v ,Pt > 0. As for the mass parameter, we take m = O.la^ 

The initial conditions are determined as follows. In the noninteracting theory, and in the 
expanding coordinates, the statistical function takes the generic form (see Appendix) 

F (r,r',p) = ^Re{H^(m T r)H^(m T r')} Q + n)j , (15) 
wherein m T = \/p\ + m 2 and 

n P 5 {2 \p - q) = (a ] p a q ) , (17) 



The bare mass m is negligible almost everywhere when solving the evolution equation. The only reason we 
keep Tji finite is to avoid an IR divergence in the integral (I2TI 1 in our renormalization prescription, J dprFo ~ 

2 In actual simulations, the transverse momenta squared is replaced by 

pi 2 / irriT \ 1 / irriT \ 

5-(^F( 1 - COS ^2 I ) = 2( 1 - COS ^)' (16) 
in order to reduce lattice artifacts. 
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is an arbitrary distribution of pre-existing quanta. We assume that the initial state contains only 
the classical field © and no on-shell particle excitations. This means n p (r ) = 0, and the 
classical field induces an effective mass term for the fluctuation. Our initial conditions are thus 
given by 

F\\{t,t', P ) = ^Re{H^{m\\r)H^{m^')) , m\\ = yj p| + m? + 3o* , 
F ± (r,r',p) = ^Re{Hff v (m ± T)H^(m ± r')Y m ± = yJp* T + m? + a* , (18) 

evaluated at r = r' = r . 

Since the equations (f8l-(fTTI) are second order in the r-derivative, we must also specify the 
first derivatives. For the classical field, we take 

«9 t 0(t)| to = O, (r ->0), (19) 

again motivated by the situation in the CGC case [Q]], whereas d T F(r, r ) | T0 and d T d T /F(r, r')\ T0 
are computed from ([TBI . The initial condition for the spectral function is fixed by the canonical 
commutation relation 

Pab(r , r , p) = , d T p ab (r, r ) | r _ rn = S ab — . (20) 

(l20l) is actually satisfied for all values of r even in the interacting theory, so it can be used as a 
check of the simulation. 

Concerning renormalization, in (l+2)-dimensions the tadpole self-energy © is linearly di- 
vergent. We renormalize this by subtracting the free part at each step of the evolution 

F(r,r)^ J ^L(F(t,t,p)-F (t,t, P )) , (21) 

where f F is computed from (fl3T ) (with n p = 0). There is also a logarithmic divergence in 
the 'sunset' diagram which appears in E at two loops. Though in principle this requires an 
extra term in the mass renormalization (see, e.g., U16|]), in this work we do not perform the 
renormalization of logarithmic divergences. 

As stated in the introduction, our main interest is the difference between the full 2PI calcula- 
tion and its classical statistical approximation. We thus show the result of the latter simulation 
as well. In the present context (at NLO in l/N), this simply amounts to neglecting, in the set of 
equations (l8l)-(fTTT). terms quadratic in p's when they appear in conjunction with terms quadratic 
in F's 1120. Note that we still perform the renormalization (|2~T1 . so in this sense our 'classical' 
simulation is not exactly the same as the other classical statistical simulations where cut-off 
dependent tadpoles are kept as they are. 

Furthermore, in order to see the effect of the expansion, we also perform simulations in the 
non-expanding case in the flat coordinates (t,xr,z) on an isotropic lattice a z = 2/ a = a^, 
a t = a T . The parameters are chosen such that the relation a z = tqcl^ (from dz = rdrj at 
midrapidity 77 = z = 0) is satisfied at the initial time t = r . The initial conditions for the 
classical field are the same as in the expanding case (up to the replacement r — > t), while those 
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for the F, p-functions are 



1 / o o 2\ 



F\\(t ,t ,p) = —======, (p = p T +p zJ , 

2yp + m 2 + 3a 2 

F±(to t ta,p) = 1 (22) 

ZyP + m + ° 

Pab{t,t,p) = , <kPabtyJ,p) \ t=t , = 8ah, (23) 

and d t F(t,t Q ,p)\t=t = 0, etc. (c.f., Q). 

The degree of thermalization may be inferred from the functional form of F at large time. 
In the non-expanding case, we extract the effective distribution n p (t) from the data using the 
prescription 

\ + n p {t) = y/F(t,t',p)d t d tl F(t,t',p) | t=t , , (24) 

which is often employed in the literature (e.g., Refs. The 'equilibrium' distribution 

in the expanding case is not known, nor is it clear if such a distribution is theoretically well- 
defined in our anisotropically expanding system. We try, however, the following ansatz for the 
late-time behavior of the correlation function (cf. Eqs. (fT3T ) and (IA.10I) ) 

C (\ 1 



F(t,t, P ) = — ( g + e^/T,! ) ' "p=yJPT + ml ff , (25) 



and see if we can identify an exponentially decaying (in p T ) component by adjusting the pa- 
rameter C(t) » 1 ('field strength renormalization'). At the same time we determine the fitting 
parameters T{j) ('temperature') and m e //(r) ('thermal effective mass'). 

Note that, as already implied by (l25l . we neglect the dependence of F on p v , or equivalently, 
the longitudinal momentum p z ~ P v / r at late times. Indeed, in our simulations F becomes 
independent of p v after some time because p v enters the equation of motion in the form (j5 r/ /r) 2 
(see, also, (IA.13I) ). In order to see the dependence on p z ~ p v /r, we have to use a considerably 
larger p ri cutoff such that the typical value of p^/r at large r is of the same order as the pr cutoff. 
We do not perform such simulations in this paper both from technical and phenomenological 
reasons. Technically, 2PI simulations in a very large (and anisotropic) lattice is computation- 
ally very expensive. Besides, the product of Hankel functions used in the initial condition (TT~8T ) 
behaves badly when p v becomes larger than around 10. [Our current cutoff is p v < 2n.] Phe- 
nomenologically, in heavy-ion collisions, experimentalists usually measure the p T dependence, 
but not the p z dependence (to our knowledge) because the latter is not boost invariant and de- 
pends sensitively on the collision energy. Our primary goal here is to reproduce the exponential 
Pt dependence around midrapidity rj ~ which is well established experimentally, and for this 
purpose the inclusion of finite-p^ modes is presumably not very crucial: Due the the redshift 
Pz ~ Pti/t, the relevant values of p z naturally decrease with time, and modes with finite p z 
escape from the midrapidity region and do not affect much the pr distribution there. 
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4. Numerical simulations: Results 



In this section we present the numerical solution of the equations (l8l)-(fTTT) obtained after 1600 
steps (t/t = 160) of evolution in time. We run 2 3 = 8 different simulations corresponding to 
the following cases: (i) weak A = 10 ~ 4 and strong A = 10 coupling, (ii) 2PI and the classical 
statistical approximation, (iii) with and without longitudinal expansion. All the results below 
are for N = 4 . 

4.1. Weak coupling 

Let us first present the results at weak coupling A = 10~ 4 . Fig.[T]shows the time evolution of 
the classical field (in units of y/a) in the expanding case (left) and non-expanding case (right). In 
the expanding case, the amplitude of decays very fast and has a characteristic shape. This can 
be understood analytically [|17|1 at weak coupling where the evolution of <p is largely determined 



by the classical equation. A detailed analysis gives 0(r) ~ 1/r ' , while the period of oscil- 
lation is constant in t 2 / 3 , meaning that the effective period becomes longer at larger r. These 
tendencies are clearly visible. In the non-expanding case, on the other hand, the amplitude is 
constant in time at the classical level, and diminishes only as a result of its coupling to quantum 
fluctuations. The difference between the full 2PI and the classical statistical approximation is 
unnoticeable, hence we did not plot the latter. 





Figure 1 : Time evolution of the classical field at weak coupling A 
expansion. 



10 4 . Left: with expansion; Right: without 



Next, in Fig|2]we plot Fn(r, t,Pt,P v = 0) (expanding case, left) and F\\(t,t,p T ,p z = 0) 
(non-expanding case, right jj for three different values of px/cr = jjgTiT , nr = 0, 8, 16. We 
confirm the parametric resonance previously observed at weak coupling [| 1 VM . We also see 
that the effect of the expansion is huge, suppressing the magnitude of F by several orders 
of magnitude compared with the non-expanding case. Actually, the suppression is so strong 
that the equal-time correlator F(r,r,p), which should be a positive quantity, turns slightly 
negative in some momentum modes beyond r > 50r . We presume this is because of large 
errors incurred when very large and very small numbers (which differ by ~ 10 6 ) coexist, or 



In the nonexpanding case, we actually plot aF(t, t,p) to make it dimensionless. 
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simply due to an artifact of the smallness of our lattice^ We note, however, that F(r,r,p) 
turns positive again for all momentum modes at later times r > 100r , and shows a smooth, 
nonthermal distribution in p?. 




Figure 2: Time evolution of the statistical function at weak coupling A = 10 for three different values of tit- 
tit = 0, 8, 16 (from top to bottom). Left: with expansion. Right: without expansion. 

Again, there is no discernible difference between the 2PI and the classical statistical approx- 
imation. This is understandable in the non-expanding case as F becomes very large and the 
condition FF 3> pp is always satisfied. [Remember that p is an 0(1) quantity due to the 
constraint (l23l).l In the expanding case, it is simply because the whole collision term is small 
and contributes very little to the evolution. In either case, the lack of any difference in the two 
simulations indicates that the system is far from the true equilibrium state even at the final time 
t/t = 160. The pj—distribution is so overpopulated in the soft region that it is entirely impos- 
sible to find any hint of 'subtle' quantum effects like the exponential distribution. Presumably 
we have to wait for a much longer time to see the difference, but this is beyond the scope of the 
present work. 

4.2. Strong coupling 

We now turn to the results at strong coupling A = 10. In Fig.[3]we plot the time evolution of 
the classical field. Compared to the weak coupling case (Fig. Q]), the decay of the amplitude is 
faster because the transfer of energy to quantum fluctuations is more efficient at strong coupling. 
Still, the difference between the 2PI and classical simulations is barely noticeable. 

Next we plot the statistical function at strong coupling in Fig. |4] which should be compared 
with the weak coupling counterpart Fig. |2] Because the magnitude of the classical field is 
smaller, the parametric resonance is less pronounced, and is basically unrecognizable in the 
expanding case. At least in the present model, this ineffectiveness of the instability appears to 
be the key to achieve an early thermalization. Indeed, we get a reasonable fit of the momentum 
distribution using the ansatz (|25l ). Fig. |5] shows the distribution n p {r) = tco p F(t, r,p)/C — | 



4 This problem does not seem to occur in the Monte Carlo simulation (rather than solving differential equations) 
of the classical statistical theory on a larger lattice lfl7ll . It is computationally very expensive and not realistic for 
us to increase the lattice size. 
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Figure 3: Time evolution of the classical field at strong coupling A = 10. Left: with expansion; Right: without 
expansion. Both the 2PI result (black solid line) and the classical statistical result (red dashed line) are shown, 
although the difference is very small. 




Figure 4: Time evolution of the statistical function at strong coupling A = 10 for three different values of tit' 
tit = 0, tit = 8, tit = 16 (from top to bottom). Left: with expansion. Right: without expansion. 

in the expanding case at two values of r : t/t = 50 and r/r = 150. The fitting parameters 
are found to be 





C 


T/a 


m eff /a 


t/t = 50 


0.905 


0.354 


0.0817 


r/r = 150 


0.913 


0.413 


0.0586 



where C ~ 1, as expected. The exponential behavior actually starts to show up already around 
r/r > 40. The difference between the 2PI and the classical statistical approximation is now 
manifest. In the latter case (red lines), the distribution plunges and becomes negative in the 
high-p T regiorjf] (around n T = 15 ~ 20) where the quantum (2PI) distribution still keeps up 
with the exponential trend. 



5 This is an expected behavior because the equilibrium distribution in the classical statistical theory corresponds 
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Figure 5: The distribution n p as a function of lo p at strong coupling A = 10 in the expanding case. Left: r = 50ro. 
Right: r = 150tq. The distribution from the classical statistical approximation (red line) turns negative beyond 



u p > 0.8. 



Similarly, in the non-expanding case we plot in Fig. |6]the distribution n p defined in (|24|) as a 
function of p = \/ p\ + p\. Again the distribution in the 2PI case develops an exponential tail 
already at t = 50t , but the classical distribution follows a power law and becomes negative. 
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Figure 6: The distribution n p defined in (|24| > at strong coupling A = 10 without expansion. Left: t — 50io; Right: 
* = 150< - 



to the following approximation to the Bose-Einstein distribution 

1 T 1 



pi 1 — 1 lj 2 
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5. Conclusions 

In this paper, we have for the first time performed the 2PI simulation of the O(N) model with 
longitudinal expansion. This model is unique in that, to NLO in the l/N expansion, one can 
resum the 2PI diagrams to all orders in the coupling constant. Therefore, it is valid for all values 
of the coupling and also in the presence of violent instabilities. At weak coupling, the parametric 
resonance is strong and the classical statistical approach is an extremely good approximation. 
This means, however, that we cannot see any hint of the quantum equilibration within the time 
range explored in this work. Just for a (naive) comparison, we note that in heavy-ion collisions 
one would choose r ~ l/Q s ~ 0.1 fm (the inverse of the saturation momentum) and the 
lifetime of the QGP is < 5 fm, so in practice we have a constraint t/t < 50. Thus, in order to 
achieve the truly equilibrated system (not just the establishment of an equation of state) in such a 
short time, pictures based on the weak coupling limit may not be a good starting point. Instead, 
we have seen that in the 2PI simulation at strong coupling the difference from the classical 
statistical approximation does show up within a reasonable time scale, and plays a crucial role 
in developing the Bose-Einstein (or the exponential) distribution in px. 
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Appendix A. Free correlation function 



The quantization of the scalar field in the t-t] coordinates was studied in [|18|1. In (1+2)- 
dimensions, the canonical commutation relation 



[(p(r,x T ,r]),Td T (f)( y T,x' T ,r]')] = i8(rj - r)')6(x T - x' T ) , 
admits the following mode expansion 

(f)(x) = j dp T dp v a PTjPv £ PTjPv (x) + al T ,p v €p T , Pv ( x ) 
Here the one-particle wavefunction 



(A.l) 



(A.2) 



f 9 

/ TdxTdrjQ^x) i-7^€ Pr ,,p T (x) = 5(p v - p' v )5(p T - p T ) , (A.3) 



with m T = \/ p\ + Tn 2 , is normalized such that the commutation relation of a and takes the 
form 



Kt,Pvi a p' T , P 'J =$(Pt- Pt)$(Pv ~ P't,) ■ 



(A.4) 
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The statistical and the spectral functions of the free theory are easily calculated 

F(r,r',p) = I(O|{0(r,p),0(r',-p)}|O> = ^e(^Kr)^(m T r')) , (A.5) 

p(r,r',p) = z(O|[0(T,p),0(r',-p)]|O) = ~Im (^(m^^r^r')) . (A.6) 
The approximate form at early times r, r' — > is 

BY , „ 1 cos(2 Pr; ln^-2%,)) 

.Mr, r, p) ps 1 =- , (A.7) 

2p r; tanh^7r 2^ sinh^vr 

p(t,t',p) ~ — sin/p^ln — ), (A.8) 

where 

0(p„) = arg(r(ip„)) , (A.9) 
is the phase of the gamma function. The late-time behavior r, r' — >■ oo is 

F(r,r',p) ~ 1 cos m r (r - r') , (A. 10) 

2mr vtt' 

p(t,t',p) ~ sin mr(r — r') . (A. 11) 



rriTVTT' 



These are similar in form to the corresponding functions in the flat coordinates 

F(t,t',p) = cosujp(t — t') , p(t,t',p) = — smujp(t — t r ), (A. 12) 

2,ujp ojp 



where u p = ^Jp 2 + m 2 . Finally we note the first subleading term at r = r' (— >■ oo) 



F(r,r,p)« — + . (A.13) 
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